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Ph Abstract 

We review a recent approach to the mean-field limits in neural networks that 

Stakes into account the stochastic nature of input current and the uncertainty 
I I in synaptic coupling. This approach was proved to be a rigorous limit of the 

network equations in a general setting, and we express here the results in a 
more customary and simpler framework. We propose a heuristic argument to 
^ derive these equations providing a more intuitive understanding of their origin. 

These equations are characterized by a strong coupling between the different 
QQ moments of the solutions. We analyse the equations, present an algorithm to 

simulate the solutions of these mean-field equations, and investigate numeri- 
' cally the equations. In particular, we build a bridge between these equations 

and Sompolinsky and collaborators approach [38l [12] , and show how the cou- 
phng between the mean and the covariance function deviates from customary 
approaches. 

X Introduction 

The problem of modeling neural activity at scales integrating the effect of thou- 
sands of neurons has been a great endeavor in neurosciences for several reasons. 
First, it corresponds to the mesoscopic scale that most imaging techniques are 
able to measure. The activity monitored by these classical imaging techniques 
(such as the electro-encephalogram, the magneto-encephalogram, the optical 
imaging, etc. . . ) shows the global behavior of a whole area that results from 
the activity of several hundreds to several hundreds of thousands of neurons. 
Second, anatomical data recorded in the cortex reveal the existence of struc- 
tures, such as the cortical columns, with a diameter of about SO/um to 1mm, 
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containing of the order of one hundred to one hundred thousand neurons be- 
longing to a few different types. These columns exhibit specific functions as a 
result of the collective behavior of several neurons, as evidenced for instance in 
the primary visual cortex VI, where cortical columns are organized in a specific 
geometric configuration and respond to preferential orientations of bar-shaped 
visual stimuli. In this case, information processing does not occur at the scale 
of individual neurons but rather emerges from the collective dynamics of many 
interacting neurons contributing to a mesoscopic or macroscopic signal. The 
description of this collective dynamics requires models which are different from 
individual neuron models. Indeed, when the number of neurons is large enough, 
averaging effects appear, and the collective dynamics is well described by an 
effective mean-field, summarizing the effect of the interactions of a neuron with 
the other neurons, and depending on a few effective parameters. 

Models including these features which use population rates or activity equa- 
tions are the foundation of a common approach to modeling for neural networks. 
These equations provide meanfield dynamics for the firing rate or activity of neu- 
rons within a network given some connectivity. This approach was introduced 
by Wilson and Cowan in the 1970s [3TJI32], and has been widely studied since 
then. This approach assumes that the activity coming from a cortical column 
is captured in the mean-firing rate of the cells in the population. 

One of the principal limitations of firing rate models is that they only take 
into account mean firing rates and do not include any information about higher 
order statistics such as correlations between firing activity of different neurons 
of the same populations or across populations. We know that spike trains of 
individual cortical neurons in vivo are very noisy and have interspike interval 
(ISI) distributions close to Poisson (see e.g. [33)- However, at the level of a neu- 
ral population, the question of the statistical distribution of the activity is still 
open, and its description has been the subject of a many studies in the last fifteen 
years. Numerous studies on integrate-and-fire networks have shown that under 
certain conditions, even though individual neurons would exhibit Poisson-like 
statistics, the neurons fire asynchronously so that the total population activity 
evolves according to a mean-field rate equation, with a characteristic activation 
or gain function [U [71 [S] ■ Formally speaking, the asynchronous state only ex- 
ists in the thermodynamic limit where the number of neurons tends to infinity 
(finite-size deviations from this limit are addressed in [8]). Moreover, even if the 
asynchronous state exists it might not be stable. This has motivated a num- 
ber of authors to analyze statistical correlations and finite-size effects in spiking 
networks [321 IMl HZl HI] • Recently, a few papers analyzed the effects of correla- 
tions by utilizing the master equation introduced by El Boustani and Destexhe 
|15j and providing different moment truncations using statistical physics tools, 
such as path integrals ^ and the Van Kampen expansion [Gj . These approaches 
enlighten a complex interplay between the mean activity and the mean correla- 
tions, that yield non-trivial qualitative effects, that begin to be uncovered (see 
[39)). A complementary approach has been developed by Cai and collabora- 
tors based on a Boltzmann-like kinetic theory of integrate-and-fire networks 
[32j ITO| , which is itself an extension of population-density methods [301 El] . A 
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recent alternative approach proposes a mathematically derived mean-field equa- 
tion on stochastic firing-rate networks [TQJ taking into account variability in the 
connectivity weights. They rigorously derive the mean-field equation of the net- 
work when the number of neurons tends to infinity and study the solutions for 
the dynamic and stationary mean-field equations. 

The present article builds upon this latter approach, and aims to relate this 
abstract approach to Wilson and Cowan deterministic mean-field equation [411 
|42] on one hand, and few other stochastic methods on the other hand. To this 
purpose, we implement a simple yet very familiar model that fits the framework 
of [in]. We propose, in this simpler case, a heuristic argument that yields 
the same mean-field equation they obtain, analyse these equations, describe an 
algorithm to compute numerically the unique solution of these equations, and 
finally, investigate numerically the solutions of these equations. 

1 Material and Method : The stochastic firing- 
rate model 

In this article we study a model that we call the stochastic firing-rate models a 
stochastic version of the customary firing rate models. This section is aimed to 
define this model. 

1.1 The firing-rate model 

The building block of the network we study here consists of firing-rate Wilson 
and Cowan models [HI |Hj , where we take into account the stochastic nature of 
synaptic inputs. In details, we consider a network composed of N neurons in- 
dexed by I S {1, ... , N} belonging to P populations indexed by a € {1, . . . , P}, 
and let us denote to be the number of neurons in population a. Each neuron 
is described by its membrane potential Vi . Alternatively, this model can be seen 
as a model of cortical area, where the interacting units are cortical columns that 
follow Wilson and Cowan equations. 

The Wilson and Cowan model supposes that each incoming spike provokes 
a postsynaptic potential when received, and that these postsynaptic potentials 
are summed linearly to produce the membrane potential Vi{t). A single action 
potential from neuron j received at time t* by the neuron i provokes the addition 
of potential PSPij{t— s), where s is the time of the spike hitting the synapse and 
t the time after the spike. We neglect the delays due to the distance travelled 
down the axon by the spikes. 

Under the assumption that the post-synaptic potentials sum linearly, the 
average membrane potential of neuron i is given by the sum of all postsynaptic 
potentials coming from the other neurons j, through the formula: 

l/,(t) =^P5P,,(t-i^^) = / P5P„(t-s)^<5(t^"-s)ds (1) 
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where the sum is taken over the arrival times tj of the spikes produced by the 
neurons j. When the number of neurons tends to infinity, or when the firing 
rate of neurons increase, the sum can be approximated by the mean firing rate 
of each neuron, to first order. More precisely, this approximation consists in 
supposing that since the number of spikes arriving from neuron j between t and 
t + dt is Vj{t)dt, we can approximate equation ([!]) by: 



.,■ Jtn 





We now assume that the neuron received input for arbitrary long times, 
i.e. = — oo. Another very important assumption of the rate model is that 
the relation between the mean firing rate and the membrane potential can be 
expressed in the form Vi{t) — Si{Vi{t)), where Si is sigmoidal function (called 
the firing rate, or squashing function, see e.g. [551 [T3], a-nd this function can be 



approximated in some cases using an averaging method (see section 1.3) 



j 

We finally note that the postsynaptic potentials PSPij can depend on several 
variables in order to account, for instance, for adaptation or learning. However, 
in this paper, we follow the common assumption, initially made by Hopfield in 
|25| , that although the sign and amplitude may vary, the post-synaptic potential 
has the same shape no matter which presynaptic population caused it. This 
leads to the relation 

PSP^J{t) - .hJg^{t). 

Qi represents the unweighted shape (called a g-shape) of the postsynaptic poten- 
tials and Jij is the strength of the postsynaptic potentials elicited by neuron j 
on neuron i. In this paper, we further assume the g-shape to be an exponential 
function, and that both the g-shape and the sigmoidal functions only depends 
on the population (a) the neuron i belongs to: gi{t) = gait) = e~*/'^°lt>o- 
Under these assumptions, the membrane potential of neuron i from population 
a can be written as the solution of the following differential equation: 

dy„ 



/* PSP^.ii 

J — OO 



This equation models the interactions between the neurons in the network. 
External inputs also have a great influence in the system. These inputs are 
modelled as the sum of a deterministic input and a noisy input described in 
section [l.2[ The deterministic input, generally assumed to originate from the 
thalamus, accounts for sensory inputs and cortico-cortical deterministic activity. 
These inputs arc assumed to be identical for all neurons in the same population, 
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and are denoted /^(t) for the input to the population a. Therefore, we end up 
with the deterministic modeh 

AY ] ^ 

which corresponds to a case of the Wilson and Cowan system. We now include 
the effect of noise in this model. 



1.2 Randomness in the network 

Cortical areas receive a constant bombardment of action potentials originating 
from different sources, ranging from sensory stimuli to motor areas and cortico- 
cortical activity. These signals are characterized by a high level of variability, 
which is a result of the randomnes^ of processes arising in neuronal phenom- 
ena (see e.g.[37| |35]). This noise can have different origins. It is the result 
of contributions attributed, for instance, to thermal noise arising from the dis- 
crete nature of electric charge carriers, the behavior of ion channels, synaptic 
transmission failures and global network effects. We also include in this noise 
term all the unrelated activity either coming from external stimuli or from ex- 
tra network cortico-cortical activity. This bombardment of "random" spikes 
results in a noisy current rii whose probability distribution only depends on the 
population the neuron belongs to. Using the classical diffusion approximation 
(see e.g. [211 [TH |40]), we model that the noisy current rn as a white noise, the 
formal differential of a centered Wiener process with standard deviation for 
neurons in population a, and we assume that the noise received by each neuron 
is independent of those received by the other neurons. 

Moreover, the synaptic weights Jij are known up to a certain precision, and 
cannot be defined individually with high precision. The mean connectivity, 
and the standard deviation of the error done in the measurement can neverthe- 
less be experimentally evaluated. The random synaptic weights are therefore 
assumed to be drawn from a Gaussian distribution of mean Jap/Np and stan- 
dard deviation aap / ^jNp, and do not evolve (they are said to be frozen during 
the evolution). The scaling chosen ensures that the local interaction process 
'^j^i JijSi3{Vj{t)), summarizing the effects of the neurons in population /3 on 
neuron i, has a mean and variance which do not depend on Njj and are only 
controlled by the phenomenological parameters Jap , (Jap ■ We finally make the 
simplifying technical assumption that the Jij are independenl|^ 

^Some authors consider this variabihty as a result of an underlying finite dimensional chaos, 
see e.g. |20l |26] 

^ It is however known that synaptic weights are indeed correlated (e.g. via synaptic plastic- 
ity mechanisms); these correlations are built by dynamics via a complex interwoven evolution 
between neurons and synapses dynamics and postulating the form of synaptic weight cor- 
relations requires, on theoretical grounds, a detailed investigation of the whole history of 
neuron-synapse dynamics. 
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Under these assumptions, the membrane potential of neuron i in population 
a satisfies the equation 



dV,{t) 



, P N„ \ 

— V^{t) + Io.{t) dt + UdW,{t) (2) 

^" ,9 = 1 , = 1 / 



1.3 Averaging Models and Firing rate functions 

The nonlinear form of the firing rate function has important implications; there 
are many choices that various authors have used. The simplest is the step 
function, where the neuron fires maximally or not at all depending on whether 
the potential is above or below threshold. If one uses a statistical-mechanical 
approach to derive 'mean-field' equations [2] then this sharp step function is 
smoothed out and can be represented by one of the well-known squashing func- 
tions: 

/^(^)= i+e-f°-^T)/v. or 

\siV) = + erf(-(l^ - Vt)/V,)) 

The first one is often referred to as the logistic function, and the second as the 
Gaussian function; Vr is an activation threshold and Vs governs the slope of the 
sigmoid. These choices appear somehow subjective and arbitrary. 

When the number of neuron increases, this squashing function can be com- 
puted explicitly in some simple cases, using averaging properties (see e.g. |34) 
for class II neurons and 16 for class I neurons). This method can guide the 



choice of the squashing function. It is described in more detail in section 2.1 



2 Theory and Calculations 

The network we consider is governed by the N equations such as ([2|. In these 
equations, the neuron i of population a integrates external deterministic and 
stochastic inputs, and receives from the network neurons a current, which we 
call the microscopic interaction process. The mean-field problem consists in 
finding an equation that governs the activity of the neurons, at the level of 
the network population, when the number of neurons tends to infinity. In the 
stochastic setting of this article, this question is understood as a limit in law, 
under the joint law of the connectivities and the Brownian motions. In this 
section, we first review a way to obtain the squashing functions, before deriving 
the mean-field equations for the networks we are interested in. 



2.1 Averaging method and squashing functions 

The squashing function corresponds to the averaged effect of many intercon- 
nections. When the number of neurons increases, this squashing function can 
be computed explicitly in some simple case, using averaging properties when 
assuming slow synapses (see e.g. [34] for class II neurons and fT6| for class 
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I neurons). This method can provide a choice of squashing function. Let us 
consider for instance the foUowing system of coupled neurons: 

<f^=q{V„W,) (3) 

The jth ceU of the network is represented by its potential, Vj, and all of the 
auxiliary channel variables that make up the dynamics for the membrane, Wj. 
The term /Jpp' is any tonic applied current. Finally, the synapses are modeled 
by simple first order dynamics and act to hyperpolarize, shunt, or depolarize the 
postsynaptic cell, e is a small parameter corresponding to the time scale of the 
synapses compared to the time scale of evolution of the membrane potential, 
and is assumed to be small. Associated with each neuron is a synaptic channel 
whose dynamics is governed by the variable Sj which depends in a (generally 
nonlinear) manner on the somatic potential. Thus, one can think of Sj as being 
the fraction of open channels due to the presynaptic potential. The functions Sj 
have maxima of 1 and minima of 0. The effective maximal conductances of the 
synapses between the cells are in the nonnegative numbers Xjk and the reversal 
potentials of each of the synapses are V^"^"^. The goal is to derive equations that 
involve only the Sj variables and thus reduce the complexity of the model while 
retaining the qualitative features of the original. 

The main idea is to exploit the smallness of s and thus invoke the averaging 
theorem on the slow synaptic equations. Each of the slow synapses is held 
constant and the membrane dynamics equations are solved for the potentials, 
Vj(t; si, . . . , s„). The potentials, of course, naturally depend on the values of 
the synapses. We will assume that the potentials are either constant or periodic 
with period T{si, . . . , s„). Once the potentials are found, one then averages the 
slow synaptic equations over one period of the membrane dynamics obtaining: 

ds ■ 

= e{Sj{si, ■■■ ,s„) - Sj) 

where 

^ .T,(si,...,s„) 

S'j(si,...,s„) = — -/ Sj{Vj(t;si,...,Sn))dt. (4) 

lj[Si, . . . , Sn) Jq 

Remark. Note that if we assume that when the cell is not firing, the potential 
is below the threshold for the synapse, we obtain a model very close to the 
activity-based (see [13 )• 

This method allows derivation of the squashing function. The main assump- 
tion is that the detailed phase and timing of individual spikes is not important; 
that is, one is free to average over many spiking events. Rinzel and Frankel 
apply the same averaging method to derive equations for a pair of mutually 
coupled neurons that was motivated by an experimental preparation. In their 
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paper, they require only cross connections with no self-self interactions, and they 
are able to numerically determine the potential as a function of the strength of 
the synapses. They use a class of membrane models that are called "class 11" 
(see |33| ) where the transition from rest to repetitive firing occurs at a subcriti- 
cal Hopf bifurcations that appears, for instance, in the Hodgkin-Huxley model. 
This latter assumption implies that the average potential exhibits hysteresis as 
a function of the synaptic drive, i.e. the functions Sj{sk) are multivalued for 
some interval of values of Sk- Because of this hysteresis, they are able to combine 
an excitatory cell and an inhibitory cell in a network and obtain oscillations, 
which is impossible for smooth Sj{sk) without self-excitatory interactions (see, 
e.g. [33]). 

Class I membranes are simpler, as proved by Ermentrout in [16 . It is known 
(see e.g. [18]) that the frequency of oscillation in Class I membranes evolves as 
the square root of the bifurcation parameter (denoted by p) : 

uj = C{p*)ViP-P*)++Oi\p-p*\). (5) 

If wc make the approximation that Sj are simple on/off function, that is equal 
to if the membrane is at rest (and thus below threshold) and 1 if the voltage is 
above threshold. Let ^{p) denote the amount of time during one cycle that the 
potential is above threshold. Then the equation Q, depending on the critical 
parameter p on which the potential depends, is simplified to 



Tip) Jo ' 2^ 

where we have used the fact that T = 2tt/uj. A fortuitous property of class 
I membranes is that the time for which the spike is above threshold is largely 
independent of the period of the spike so that £^[p) is essentially constant. (This 
is certainly true near threshold and we have found it to be empirically valid in 
a wide parameter range.) Thus, combining this with ([S]), we obtain the very 
simple squashing function: 



S{p) = C{p*)^{p-p*)+. 

where C (p) and p depend on the other parameters in the model. 

The value p* and the constant C{p*) need to be evaluated. Ermentrout in 
|16j provides closed-form expressions for these quantities in the case of Morris- 
Lecar membrane model, and obtains a good agreement with the full system in 
numerical simulations. The resulting equations have the very nice compact form 
in a two cells network: 



T, 

Ti ^ + Si = Ci \J {XeiSe ~ ffg [XaSi, li)) 



where g* (g, I) is the two-parameter surface of critical values of the excitatory 
conductance. Note that these squashing functions are not similar to the usual 



8 



nonlinear transforms between the voltage and the firing rate. They are less 
smooth than the usual ones, and in particular, are not differentiable at p = p*: 
the slope at the bifurcation is infinite. The differential equations, such as ([6]), 
obtained with this squashing function do not satisfy the Lipschitz standard con- 
dition for the existence and uniqueness of solutions. Moreover, in one dimen- 
sion, differential equations with a square-root nonlinearity are classical counter- 
examples to the uniqueness condition: there can exist multiple solutions starting 
from the same initial condition. Therefore, this mean-field method yields equa- 
tions that are complex to investigate, both in the class I case where the function 
is multivalued as in the class II case where it is non-smooth. Note eventually 
that this approach yields a mean-field description of model ([s]) in the sense that 
it results from an averaging of the behavior of an assembly of cells, when the 
number of neurons tends to infinity. 



2.2 Derivation of the Mean-Field equations 

Considering now a network of firing rate neurons of type (H, a rigorous proof 
of the existence of the mean-field limit and the characterization of this limit 
is provided in J19], where the authors derived their calculations from the work 
of Ben-Arous and Guionnet [S]. This proof is extremely technical, but can 
be approached with non-rigorous heuristic arguments. We develop the heuristic 
proof of the derivation of this mean-field equation because we believe it provides 
a better understanding of what the technical mean-field limit derived in |19) 
physically corresponds to. 

The heuristic proof we propose here is based on Amari's local chaos hypoth- 
esis introduced in 1972 and then widely used |2 El ES [HI |TTJ and can be 
expressed as follows: 

Amari's Local Chaos Hypothesis: For N is sufficiently large, all the Vi are 

pairwise stochastically independent, are independent of the connectivity param- 
and have a common distribution population per population. 

Proposition 2.1. Under Amari's local chaos hypothesis, the microscopic in- 
teraction term U^^ corresponding to the current from population /3 on a neuron 
i in population a defined by: 

C/5:=^J.,^^(y,(i)) (7) 
i=i 

has a distribution that only depends on the populations a and /3, and converges 
in law (under the joint probability of the connectivity weights and the voltages) 
when the number of neurons tends to infinity towards an effective interaction 
process C/Jf/j such that: 

(e [U^pit)] = J^pE[Sp{Vpm 

\cov(C/J^(t),(7^^,(s)) = al^E[Sp{Vp{t))Spms))]K=^,0=s ' 



9 



where Va corresponds to the common law of the neurons in population a ensured 
by the local chaos hypothesis, and 1a is the indicator function of the set A. 

For the sake of simplicity, for any P-dimensional stochastic process X, we 
introduce the notation m^{t) = ^[Sp{Xp{t))\ and 

'sp{Xp{t))Sp{XfM) 

We now turn to the heuristic proof of this proposition. We consider that N is 
large enough for the local chaos hypothesis to be valid. Under this assumption, 
the voltages are independent, identically distributed and independent of the 
individual connectivity weights Jy . 

The microscopic interaction process ([T]) is therefore the sum of independent 
identically distributed random variables. The functional central limit theorem 
applies to our regular case, following [3T] (very weak regularity conditions are 
necessary on the summed processes to get this convergence), provided the con- 
vergence of the two first moments of the sum. 

The problem therefore reduces to the computation of the limits of the mean 
and standard deviation of the microscopic interaction process when the number 
of neurons tends to infinity. First of all, the mean of the microscopic interaction 
process for TV large enough is constant, since we have 

E(c/^^(t)) =E(^J,,5^(T/,(t))) =5^E(j,,)E(5^(y,(t))) 

This constant value is thus the limit of the means of the sum of independent 
random variables, and therefore the mean of the limiting Gaussian process. 

Similar straightforward computations (see appendix |A]) prove that the co- 
variance of the microscopic interaction process converges to the expression given 
in the proposition. From the functional central limit theorem and the conver- 
gence of the two first moments of the microscopic interaction, we conclude to 
the convergence of this sequence of processes to the effective Gaussian process 
of mean and covariance given in the proposition. The statistical independence 
between this effective interaction process and the membrane potential V is also 
proved in the same appendix \K\ 

As a direct consequence of this last proposition, the mean-field equation for 
the membrane potential processes can be derived. 

Proposition 2.2. Under the local chaos hypothesis, the network equations ([2| 
converge in law towards the solution of the non-Markovian stochastic equation: 

(0 = I - — 14 (<) + + (^) I + /"'^^t" • 

V J 

where the processes {Wa{t))t>to are independent Wiener processes and [t) = 
{U^p{t); a,/? € {1, . . . , P})t is the effective interaction process. 
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Proof. The solution of the network equation ([2| with initial condition V{0) at 
t = can be written as: 

Because of the convergence in law of the microscopic interaction process to 
the effective interaction process, we have the convergence in law of the inte- 
gral term e^*~*'/'^° J7^(s) ds towards the effective term e^*"*'/"^" [/^(s) ds, 
and therefore, for any neuron i of population a, the potential converges in law 
towards the solution Va of the stochastic fixed-point equation: 

V^it) = T4(0)e-*/^- + E /* e(^-*)/^-[/,^^(s) ds + f e'^-*)/^" dW^ (10) 
p^i^o Jo 

which is equivalent to equation ([9]). □ 

The two equivalent equations Q and ( 10 ) are referred to as the mean- field 
equations. In the limit where the number of neurons tends to infinity, all the 
neurons have the same distributions, behave independently and for any neuron 
in population a, its membrane potential is solution of this set of equations. 

Remark. All the derivations have been done based on Amari's local chaos 
hypothesis, which is a non-rigorous assumption. However, the exact same result 
is obtained by rigorous mathematical methods [HI IH |24] . 

We discuss the form of these mean-field equations and the existence and 
uniqueness of their solutions in the following section. 



2.3 Analysis of the Mean-Field equations 

The mean-field equations derived in the previous section are rather unusual 
equations. Indeed, they involve the effective interaction process term, which is a 
functional of the solution of the equation. More precisely, it is a Gaussian process 
whose law depends on the statistics of the solution of the equation in a non 
Markovian way. Therefore, this equation cannot be considered as a stochastic 
differential equation. It is an equation set in the space of stochastic processes, 
or equivalently an equation on the probability distribution of the mean-field 
solution. Similar issues were observed in the case of spin glasses The 
existence and uniqueness of solutions of the dynamic mean-field equation (i.e. 
starting from a Gaussian initial condition) is proved in [19j . by considering the 
mean-field equation as a fixed point equation in the set of stochastic processes. 
Conditions are also given for a stationary solution to exist. 

The approach is based on the fact that solutions having a Gaussian ini- 
tial condition V{0) (or for stationary solutions) are necessarily Gaussian pro- 
cesses, and therefore characterizing the mean fi^ (t) and the correlation function 
C^^(s,t) is sufficient to identify the solution of the mean-field equations. How 
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do we compute these two deterministic functions? Let us start by writing the 
equations they satisfy. By straightforward computations on equations ^ and 
(10), we obtain the equations for the mean ii^ {t) and the correlation Cap{s,t) 
for s < t: 



dt 



CLit, s) = e-(*+^)/-° [var(K.(0)) + ^ (e^ - l) 

+ E^=i <p lo lo Ar,(«, v)dudv , 

[Co,p{s,t) =0 for 



11 



where Dx — exp{"X /2) /\/2tt dx is the standard Gaussian measure and 



I ^ / „ . ^ y—f== + f^p ( 



X Sp {y^Cf^{v,v)+fi^{v)^ DxDy 



(12) 



In the more general setting of flP, it is proved that these equations have a 
unique solution that can be computed via the iteration of a map T defined on 
the set of continuous mean and covariance functions. This function transforms 
a Gaussian process X with mean /i^ and covariance into the Gaussian 
process Y = -FiX) with mean /x^ and covariance given by: 



t P 

(t-s)/r. 



MrW-A^f(0)e-*/^°+ / e 
Jo 



(^ J„^E [S0{Xp(s))] + Io,{s))ds 

13=1 



13=1 



+ 00 



'p (x^vf{s) + y.^{s)) Dxds. (13) 



where we denoted v-^is) the standard deviation of X^ at time s, instead of 
Caa{s,s), and 



CL(t, .) = e-(*+^)/- [.f (0) + ^ (e^ - l) 



p 

E 

/3=1 



^2 



Jo 



e("+")/"°Af.(u,w)dMdJ, (14) 



where A;^o is defined in equation (12). 
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These equations take a particularly simple form in the case where the firing- 
rate sigmoidal functions are erf functions, where erf (y) = Dx is the repar- 
tition function of the standard Gaussian distribution. Let us consider that 
Sa{x) = erf ((7Q,a; -I- 7^) where erf is the repartition function of the Gaussian. In 
that case (see appendix [b] for the calculations) , the mean-field equations take 
the form: 

and the equation governing the standard deviation of the process reads: 



Va{t) 



■^0(^0) 



2 ^ 



P 

E 



2 



^(u+u)/t 



vp(t)+g^{vp{t)vp{s)~Cpp{t,sY) 



Dy. 



(15) 



These simpler equations allow for more efficient numerical simulations, and 
an easier analytical treatment, and will be discussed in the Results section. 



2.4 Numerical simulations algoritm 

The mean-field equations can be seen as a fixed point equation in the space 
of stochastic processes, or equivalently as a fixed point equation on the mean 
and the covariance function. In (TH], the authors show that this solution can be 
obtained through the iteration of a function that operates on the mean and on 
the covariance function. This is the approach we chose to numerically compute 
the solutions of the mean-field equations. 



2.4.1 Description of the algorithm 

In details, let X be a P-dimensional Gaussian process of mean /i"^ = {l^-a {t))a=i...p 
and covariance C'^ — {C^p{s,t))a^pe{i...P}- We use the map which makes 
the mean-field equation a fixed-point equation, as described in section |2.3[ 



From equation ( 13 ) , we can see that knowing (s) , s £ [0, t] we can compute 
Ma (t) using a standard discretization scheme of the integral, with a small time 
step compared with Tq and the characteristic time of variation of the input 
current Alternatively, since /i^ satisfies the differential equation: 



dt 



/3=i 



we can compute faster and with a better accuracy the solution using a Runge- 
Kutta algorithm. 



The covariance variable of the image of T (from equation (14)) can be 
split into the sum of two terms: the external noise contribution C^^it^s) = 



13 



-(t + s)/T 



2 



, where OU stands for Ornstein-Uhlcnbeck, 
and the interaction between the neurons. The external noise contribution is a 
simple function and can be computed straightforwardly. To compute the interac- 
tions contribution to the standard deviation we have to compute the symmetric 
function of two-variables: 



e("+^)/^''A^Ju,v)dudv, 



from which one obtains the standard deviation using the formula 

p 

To compute the function H^p{t,s), we start from t ~ and s — 0, where 
iJ^(0, 0) — 0. We only compute H^i^{t, s) for i > s because of the symmetry. 
It is straightforward to see that: 



HL{t + dt,s)^H^Jt,s) 



dt 



D^p{t,s)dt + o{dt), 



with 



— s It 



^/^-A^Jt,v)dv. 



Hence computing H^p{t+dt, s) knowing H^p{t, s) amounts to computing Dafi{t^ s). 
Fix i > 0. We have Dap{t, 0) = and 

ds 

D^p{t, s + ds)^ D^p{t, - — ) + A^p{t, s)ds + o{ds). 

This algorithm enables us to compute H^i^(t, s) for t > s. We deduce H^p{t, s) 
for t < s using the symmetry of this function. Finally, to get the values of 
H-^p{t^ s) for t = s, we use the symmetry property of this function and get: 

2dt' 



H^M + dt,t + dt)^H^Jt,t) 



+ 2D^f^{t,t)dt + o{dt). 



2.4.2 Analysis of the algorithm 

Convergence rate Using the fact that the iterations of the algorithms pro- 
duce a Cauchy sequence and using evaluations done in the proofs of existence 
and uniqueness of solutions provided in [19j , one can show that the precision of 
the algorithm is given by 

||C-^"W_pXM^||^ <^2(iV + T-to)dt + i?Ar(^) (16) 

where J^^(X) is the A^th iterate of the map T on the initial process X and 
Xmf the unique fixed point of i.e. the mean-field solution. In equation 



(16), 6*2 denotes a constant depending on the parameters of the model, Rn{x) 
is the exponential remainder, i.e. Rn{x) = X^^at a;"/n! and k is some constant 
depending on the parameters of the model. 
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Complexity The complexity of the algorithm depends on the complexity of 
the computations of the integrals. The algorithm described hence has the com- 
plexity 0{Niter{-^)^) where iViter is the iteration number of the algorithm. 



3 Results 

In this section, we compare the results provided by the present approach to 
a variety of deterministic and stochastic approches. We will be particularly 
interested in comparing the solutions of the present mean-field equation to the 
deterministic Wilson and Cowan equations [HI |42] : 



1 P 

-V^it) + J2Jc.pSpiVpit)) 



dt 

13=1 



If Ci3i3{t,t) = for all i e R, then the equation the mean-field equations 



( 11 ) reduces to: 



P 

dt Ta ^ ^ 

13=1 



Ja^Sfj ifJ-l3(t)) + Ia(t), 



which is precisely the customary Wilson and Cowan equations corresponding to 
the system where no randomness is considered neither in the input nor in the 
coefficients. This is precisely the limit for /„ = and ct^^ = for all a and /3 
of the mean-field equations. Therefore, the Wilson and Cowan equations can 
be seen as the deterministic limit of the general mean-field equations presented 
in this paper. Moreover the solutions of the general mean-field equations in the 
small noise case converge to the solutions of these deterministic equations. This 
convergence is to be taken in a weak sense, since the conditions for existence 
and uniqueness of the solution of the mean-field equations are not satisfied in 
the case where no noise is present (the setting for existence and uniqueness 
necessitates to stay in the space of continuous mean and covariance functions). 

However, when noise is considered, the mean of the solutions of the mean- 
field equations differs from the solutions of the Wilson and Cowan's system, and 
non-trivial differences appear. The precise study of the differences between the 
mean-field system and customary approaches is still an active ongoing research 
field, but evidences of non-trivial effects are already found. First of all, in the 
one-dimensional case, we obtained numerically that presence of noise can delay 
the apparition of bifurcations, stabilizing fixed points that were unstable in 
the corresponding Wilson and Cowan equation (see figure [T]). This numerical 
observation is under analytical investigation. Other non-trivial effects include 
transitions between a deterministic and a chaotic behavior, similar to the ones 
predicted analytically by Crisanti, Sompolinsky and colleagues [351 112] j and 
that we study in section [O] Another non-trivial effect is the influence of the 
correlations in the oscillatory regime for a two populations case, presented in 
[321 
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Figure 1: Delayed pitchfork bifurcation numerically simulated in a one popula- 
tion network. We plotted the value of the mean of the process on an interval 
[Ti,T2] with Ti sufhciently large, as a function of the slope g of the sigmoid 
Six)^l/{l + exp{-gx)) 

3.1 Stationary solutions in the one population case 

The mean-field equations ([9]) provide another approach to the analytical re- 
sults of the extensive work of Crisanti, Sompolinsky and colleagues [551 US]- In 
these papers, the authors use a generating functional approach to derive the 
mean- field equations the system satisfies. The equations they obtain are very 
similar to ours, in the one population case. Hence so far we have exactly the 
same equations and formalism. They further simplify the problem by consid- 
ering stationary solutions. In this particular case, the necessary conditions to 
get a stationary solutions to the mean-field equations stated in [TS] are satis- 
fied, namely, since the time constants of evolution for each population and the 
standard deviation of the noise are constant, and what they call the leak matrix 
L, here the diagonal matrix with diagonal element — I/tq, has full rank and 
strictly negative eigenvalues. The stationary solutions can therefore be written 
as: 




In the stationary case, the mean is constant, and the covariance Caait,s) only 
depends on the time difference t — \t — s\. They then write the fixed point 
equation on the covariance as a second order differential equation. This equation 
can be written as the mouvement of a particle in a potential well, this potential 
depends on a non-free parameter, q, which is equal to the standard deviation of 
the process (for a time difference t = 0), namely C(0). Because of the structure 
of stationary solutions, this parameter depends on the whole past of the solution 
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Figure 2: Simulations of the correlation equation in the stationary case. ti = 
0.25, (Ti.i = 1, S{x) = tanh{gx), and J = 1. (A) Deterministic case: g = 0.5. 
The graph shows the correlation for 1 (red), 2 (green), 3 (blue) and 4 (purple) 
iterations of the map in a semi-log scale, and illustrate the convergence of 
the iterations towards zero as predicted by Sompolinsky et al.. (B) Chaotic 
case: g = 5. We observe after convergence of the algorithm that the covariance 
function is non- vanishing (red crosses), and is well approximated by a fourth- 
order Taylor expansion of the theoretical solution (blue dotted curve) and even 
better by the sixth order expansion (purple curve). 



on (— oo, 0), and the equation they write, of type: 

-d^-~w ^^^^ 

where Vq is the potential depending on q — C(0) , remains a functional equation 
on the covariance, that do not differ from ours. Note that this dependency 
illustrates the non-Markovian nature of the problem. 

At this point, the authors treat the problem as if q were a free parameter, 
which allows them to study the phase portrait of these equations and distinguish 
two regimes depending on the maximal slope g of the sigmoidal firing function: 
for small values of 17, a deterministic regime where the covariance has a unique 
fixed point equal to 0, and for larger values of g, a homoclinic trajectory con- 
necting the point q = C* > Q where Vq vanishes to the point C — Q. The 
authors show that this is the only stable solution in the system, and interpret 
this solution as a chaotic solution in the neural network. Moreover, this solu- 
tion can be found using energy conservation, and Taylor expansion of Vq can 
be provided. These two predicted solutions are found by the simulation of the 
full mean-field system, as illustrated in Figure [2] In this figure, we represent 
the correlation function C(s, t) as a function oi t — t — s. Hence for each value 
of T correspond several correlation values which are superposed if the solution 
is stationary (resulting in a curve). These numerical results are compared to 
Sompolinsky's predictions and to the Taylor expansion to the fourth and sixth 
order. The authors of |38t JJj predict that there exist a critical value of the 
slope of the sigmoid function 5, denoted gc^ at which the system presents a 
sharp transition. This value in the numerical case we treat in figure [2] is gc = 4. 
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For any g smaller than the critical value gc, the system will present a deter- 
ministic behavior, with a covariance function identically equal to 0. They call 
this solution the trivial case, and in that case the solutions of the mean-field 
equations reduce to the solutions of Wilson and Cowan system. Though they 
do not prove that the system indeed converges to this solution, they show that 
it exists when the parameter q is free, this solution exists. We numerically com- 
pute, using the algorithm presented and integrating from a very large negative 
time, that indeed this solution is the one that is selected by the system, i.e. 
the only solution to the stationary mean-field equation. Similarly, for g > gc-, 
the covariance of the system converges towards the convergence corresponding 



to the homoclinic orbit of equation ( 18 1 corresponding to the case where the 



parameter q = C{0) is free. Therefore, the behaviors predicted by the studies of 
Crisanti, Sompolinsky and colleagues are confirmed by our numerical approach, 
and their approximation is therefore a posteriori validated. 

3.2 Oscillations in a two-populations network 

Let us now present a case where the fluctuations of the Gaussian field act on 
the dynamics of ^a{i) in a non trivial way, with a behavior departing from 
the naive mean-field picture. We consider two interacting populations where 
the connectivity weights are Gaussian random variables Ja^ = M{Jap,<7) for 
(a,/3) S {1,2}^. We set Sij{x) = tanh(ga;) and = 0, /q, = 0, Tq, = r, a = 
1,2. Note that in that case, the nonlinear function S{x) takes negative values. 
This case has a theoretical interest in the sense that analytical calculations can 
be pursued, which provides a good understanding of the mechanisms. Note 
that this case corresponds to a case where the sigmoidal transform is positive, 
by considering that S{x) — tanh^gx) + 1 (which is always positive), and for 
which the current la is now no more constant but equal to the random variable 
—Jap- All the computations can be done similarly when the current is a random 
variable. The dynamic mean-field equation for ^a{i) is given, in differential 
form, by: 



dt Tn 



^Ja/} S (^vp{t)x + np{t)\ Dx, a = 1,2. 



/3=1 

Let us denote by Gq(/x, C{t)) the function in the righthand side of the equality. 
Since S is odd, S{^/^J^x)Dx = 0. Therefore, we have Ga(0,C(i)) = 
whatever C(t), and hence the point fii = 0,fi2 =0 is always a fixed point of 
this equation. Let us study the stability of this fixed point. To this purpose, we 
compute the partial derivatives of Ga(/i, C(t)) with respect to /i^ for {a, (3) G 
{1,2}2. We have: 

|^(/x,C(i)) = -^+5Ja/j 1" (^l-tanh^ (^,f^)x + fip{t)^'^ Dx, 
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and hence at the point fj.i = 0, 1-P2 = 0, these derivatives read: 

where h{vp{t)) = 1 — ta,TL\i'{^Jv0{t)x)Dx. 

In the case Uq(0) = 0,f7 = 0, /q =0, implying UQ,(t) = for t > 0, and the 
equation for reduces to: 



dt Tr. 



2 

-^Jaf}S{}ip{t)) 

13=1 



which is the standard Wilson and Cowan system where Gaussian fluctuations 
are neglected. In this case the stability of the fixed point = is given by the 
sign of the largest eigenvalue of the Jacobian matrix of the system that reads: 



T2 



r' \ J-. 
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^22 



The eigenvalues are in this case — ^ + gXi^2, where Ai^2 are the eigenvalues of 
S and have the form: 

, _ J11 + J22 ± - ^22)2 +4J12J21 

Al 9 — . 

2 

Hence, they are complex whenever J12 J21 < "("^ii ~ ^22)^/4, corresponding 
to a negative feedback loop between population 1 and 2. In that case they have 
a nonzero real part only if Jn + J22 ^ (self interaction). This opens up the 
possibility to have an instability of the fixed point (/i = 0) leading to a regime 
where the average value of the membrane potential oscillates. This occurs if 
•^11 + •^22 > and if g is larger than: 



9c 



t{Jii + J22) ■ 



The corresponding bifurcation is a Hopf bifurcation. 

The presence of this bifurcation and these oscillations depends on the pres- 
ence of fiuctuations of the Gaussian field. These oscillations can simply disap- 
pear, as we study in depth in a forthcoming paper. Indeed, similarly to the case 
of the pitchfork bifurcation (Figure [T]) , the covariance of the noise will globally 
have the effect to delay, or event destroy the Hopf bifurcation, stabilizing the 
fixed point /i = 0. 

We show here that in a case where the cycle is conserved, the correlation will 
interact with the mean in an intricate fashion. We observe that the correlations 
present oscillations four times faster than the mean. The mean asymptotically 
oscillates at the same frequency as Wilson and Cowan system, and in opposition 
of phase. 
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Figure 3: Evolution of the mean and variance vi{t) for the mean-field of 

population 1, for J = and J = 2, over a time window [0, 20]. n is the number 
of iterations the map J^. This corresponds to a number of iterations for which 
the method has essentially converged (up to some precision). Note that vi{t) 
has been magnified by a factor of 100. Though Gaussian fluctuations are small, 
they have a strong influence on i^i{t). 



4 Conclusion 

We reviewed a recent approach to the mean-fleld limits in neural networks that 
takes into account the stochastic nature of input currents at each level as well 
as the uncertainty in synaptic coupling. One of the main advantages of this 
approach is that it can be rigorously derived from first principles using advanced 
mathematical tools in the field of probability, such as large deviation techniques 
and fixed points in the set of stochastic processes. However, in the simpler case 
treated here, the solutions are Gaussian and therefore described by the mean 
and covariance functions. Equations on the mean and the standard deviation 
are written in the simplest possible way, and a simulation algorithm is provided. 

These equations are characterized by a strong interplay between the mean 
and the covariance functions, that makes the solutions deviate from the cus- 
tomary Wilson and Cowan system. This coupling between the mean and the 
covariance was recently investigated independently in mean-field limits of spik- 
ing neurons using a Markovian formalism, but the equations involved appear to 
be simpler, and described by ordinary differential equations [HIE]. The more 
complex integral equations that we obtain here originate from the complex mem- 
ory effects that arise in mean-field limits; this has already been observed and 
discussed in other fields, for instance in the case of spin glasses and in the 
queueing theory. This phenomenon explains why the equations obtained are 
not Markovian. The links between the present mean-field equations and the 
mean-field equations studied by Bressloff, Buice, Cowan and colleagues is cur- 
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rently an important area of investigation. Another perspective is the precise 
study of mean-field equations and of how the results depart from Wilson and 
Cowan's system. In particular, calculations can be developed in a few par- 
ticular cases, and this phenomenon provides a first approach to bifurcations 
in stochastic systems, which is a great endeavor in mathematics, physics and 
biology today. 



A Convergence of the moments of the micro- 
scopic interaction process 



In the proof of proposition |2.1[ we omitted for the sake of compactness the 
calculations necessary to identify the covariance of the effective interaction pro- 
cess and its asymptotic independence with respect to the voltage process. The 
simple yet heavy computations are provided here. 

We start by computing the covariance between the random variables [/^ (t) , U^^g 
for s > and (a, /3, 7, (5) e {!,..., P}'^, assuming that TV is big enough so that 
Amari's local chaos hypothesis applies. In that case, we have: 

Cov(t/„^^(0,f/,^,(s)) 

This product involves four types of terms. However, because of the independence 
of all involved processes, this sum simply reads: 




^a— 7,/5— (5 



We therefore proved the existence of a limit to the covariance of the pro- 
cesses, which proves the convergence of the microscopic interaction process to 
the effective interaction process, by application of the functional central limit 
theorem. 

In order to prove that the effective interaction process is independent of the 
membrane potential V, we use the fact that both processes are Gaussian, and 
therefore that independence is equivalent to a null correlation between the two 
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processes. Let us consider a neuron i of class 7, a,/3 G {1,...,P}, and k a 
neuron of population a. We have: 



= la=^^Cov(T4,^(K)). 

which tends to zero when the number of neurons tends to infinity since the 
covariance can be easily proved to be bounded. 

B Case of erf firing-rate transformations 

In this appendix we provide the computations leading to simplify the mean-field 
equations when the sigmoidal functions are erf functions Si{x) = erf(gp(j)a; + 



7p(j)). In that case, the mean equation (13 1 



27r 

dxdx' 



/2 

dx 



Let us write this equation with generic boundaries 

-dxdy 



ax+b ^^{x^+y^)/2 



R J -CXI 



2tt 



The integration domain has therefore an afhne shape as plotted in figure |4j We 
make the change of variables consisting of a rotation of the axis and making 
the boundary of the domain parallel to one of the new axis (rotating from 
{x,y) to {u,v), see figure |4]). The rotation is therefore of angle a such that 
tan{a) = ga. The new integration domain is in the new coordinates given by 

^ ^ bsin(a) _ gb 



e 



8 b 



dxdy ^ I I -^^^ e''^^"+-'^/^^dudv 



'RJ-00 27r jRj-00 27r 

gb \ 



erf 



^/lT~g^ 
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Figure 4: Change of variable for the erf function. 



which reads with the parameters of the model: 

rn^W-erff^^^');f)+;f ) (19) 

The correlation function is simplified along the same lines. Indeed, the term 
A?(t, s) can be written as: 



/ S{a X 


+ by + 




/ S{dy- 




IR 


J'R 


[ S{dy~ 





by + c 

which, in terms of the model's parameters, read: 



Dy 



7r V yy^(^t)+g^{v0{t)v0{s)-C^p{t,sy) 



jDy 

(20) 
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